Single-Cell Omics: Assignment 2

Part 3: Cell-cell Communication using CellChat

Muhammad Ali (ma8308), Shyanon Rai (rais03), Lucia Galassi (lg3680)

---

Setting up Python Environment for CellChat

library(reticulate)
use_python("/opt/anaconda3/bin/python", required = TRUE)  # Force Anaconda Python
#reticulate::py_install("umap-learn", pip = TRUE)
reticulate::py_module_available("umap")  # Should return TRUE

## [1] TRUE

reticulate::py_config()

## python:         /opt/anaconda3/bin/python
## libpython:      /opt/anaconda3/lib/libpython3.12.dylib
## pythonhome:     /opt/anaconda3:/opt/anaconda3
## version:        3.12.7 | packaged by Anaconda, Inc. | (main, Oct  4 2024, 08:22:19) [Clang 14.0.6 ]
## numpy:          /opt/anaconda3/lib/python3.12/site-packages/numpy
## numpy_version:  1.26.4
## umap:           /opt/anaconda3/lib/python3.12/site-packages/umap
## 
## NOTE: Python version was forced by use_python() function

Loading the libraries

suppressPackageStartupMessages({
library(tidyverse)
library(data.table)
#library(Matrix)
#library(hdf5r)
library(lobstr)
#library(R.utils)
library(parallel)
library(irlba)
library(ggplot2)
library(uwot)
library(patchwork)
library(Seurat)
#library(SeuratData)
#library(Azimuth)
#library(SingleR)
library(scIntegrationMetrics)
library(harmony)
library(ggsci)
library(destiny)
library(SeuratWrappers)
library(scIntegrationMetrics)
library(reshape2)
library(CellChat)
library(ggalluvial)
library(NMF)
library(umap)
})

3.1. Creating a CellChat Object

We start by creating a CellChat object from our Seurat object.

merged_seurat = readRDS("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/merged_seurat.rds")

merged_seurat@meta.data$Malignancy_Cell_Types = merged_seurat@meta.data$Coarse_Cell_Type

merged_seurat@meta.data$Malignancy_Cell_Types[merged_seurat@meta.data$PredictionRefined == "malignant"] = "Malignant"

head(merged_seurat@meta.data[,c("PredictionRefined", "Coarse_Cell_Type", "Malignancy_Cell_Types")], 10)

##                                    PredictionRefined Coarse_Cell_Type
## AML1012-D0_AML1012-D0_AAAAAGTTACGT            normal              GMP
## AML1012-D0_AML1012-D0_AAAACACCAATC         malignant              GMP
## AML1012-D0_AML1012-D0_AAAATAGCCTTT         malignant             Prog
## AML1012-D0_AML1012-D0_AAACATTAAACG         malignant          ProMono
## AML1012-D0_AML1012-D0_AAACCACGTGCN         malignant             Prog
## AML1012-D0_AML1012-D0_AAACGCTGGAAN         malignant              GMP
## AML1012-D0_AML1012-D0_AAACTGGTACTA            normal              GMP
## AML1012-D0_AML1012-D0_AAACTTGCCCGT         malignant              HSC
## AML1012-D0_AML1012-D0_AAATGAAAGTCC         malignant              GMP
## AML1012-D0_AML1012-D0_AAATGTACGGTA            normal              GMP
##                                    Malignancy_Cell_Types
## AML1012-D0_AML1012-D0_AAAAAGTTACGT                   GMP
## AML1012-D0_AML1012-D0_AAAACACCAATC             Malignant
## AML1012-D0_AML1012-D0_AAAATAGCCTTT             Malignant
## AML1012-D0_AML1012-D0_AAACATTAAACG             Malignant
## AML1012-D0_AML1012-D0_AAACCACGTGCN             Malignant
## AML1012-D0_AML1012-D0_AAACGCTGGAAN             Malignant
## AML1012-D0_AML1012-D0_AAACTGGTACTA                   GMP
## AML1012-D0_AML1012-D0_AAACTTGCCCGT             Malignant
## AML1012-D0_AML1012-D0_AAATGAAAGTCC             Malignant
## AML1012-D0_AML1012-D0_AAATGTACGGTA                   GMP

data.input = merged_seurat[["SCT"]]@data # normalized data matrix
# For Seurat version >= “5.0.0”, get the normalized data via `seurat_object[["RNA"]]$data`
labels = Idents(merged_seurat)
meta = data.frame(labels = labels, row.names = names(labels)) # create a dataframe of the cell labels

cellchat = createCellChat(object = merged_seurat, 
                          group.by = "Malignancy_Cell_Types",
                          assay = "SCT")

## [1] "Create a CellChat object from a Seurat object"
## The `meta.data` slot in the Seurat object is used as cell meta information

## Warning in createCellChat(object = merged_seurat, group.by = "Malignancy_Cell_Types", : The 'meta' data does not have a column named `samples`. We now add this column and all cells are assumed to belong to `sample1`!

## Set cell identities for the new CellChat object 
## The cell groups used for CellChat analysis are  B, cDC, CTL, earlyEry, GMP, HSC, lateEry, Malignant, Mono, NK, pDC, Plasma, ProB, Prog, ProMono, T

3.2. Setting up the CellChat Database

Before we can use the CellChat object, we need to set the ligand-receptor interaction database and identify over-expressed ligands and receptors.

The database CellChatDB is a manually curated database of literature-supported ligand-receptor interactions in both human and mouse. Some of the highlights are as follows:

• ~3,300 validated molecular interactions
• ~40% of secrete autocrine/paracrine signaling interactions
• ~17% of extracellular matrix (ECM)-receptor interactions
• ~13% of cell-cell contact interactions
• ~30% non-protein signaling

It should be noted that for molecules that are not directly related to genes measured in scRNA-seq, CellChat v2 estimates the expression of ligands and receptors using those molecules’ key mediators or enzymes for potential communication mediated by non-proteins.

When analyzing human samples, use the database CellChatDB.human; when analyzing mouse samples, use the database CellChatDB.mouse.

CellChatDB categorizes ligand-receptor pairs into different types: - “Secreted Signaling” - “ECM-Receptor” - “Cell-Cell Contact” - “Non-protein Signaling” (Not used by Default)

CellChatDB = CellChatDB.human # use CellChatDB.mouse if running on mouse data

showDatabaseCategory(CellChatDB)

The function subsetDB allows us to filter the CellChatDB. And when you use subsetDB without any parameters, it gives you all the database except non-protein signaling (which is metabolic and synaptic).

# use a subset of CellChatDB for cell-cell communication analysis
#CellChatDB.use = subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation") # use Secreted Signaling

# Only uses the Secreted Signaling from CellChatDB v1
#  CellChatDB.use <- subsetDB(CellChatDB, search = list(c("Secreted Signaling"), c("CellChatDB v1")), key = c("annotation", "version"))

# use all CellChatDB except for "Non-protein Signaling" for cell-cell communication analysis
# CellChatDB.use <- subsetDB(CellChatDB)


# use all CellChatDB for cell-cell communication analysis
# CellChatDB.use <- CellChatDB # simply use the default CellChatDB. We do not suggest to use it in this way because CellChatDB v2 includes "Non-protein Signaling" (i.e., metabolic and synaptic signaling). 

# set the used database in the object
cellchat@DB = subsetDB(CellChatDB)

3.3. Identifying Over-expressed Ligands and Receptors

cellchat = subsetData(cellchat) # This step is necessary even if using the whole database

cellchat = identifyOverExpressedGenes(cellchat)
cellchat = identifyOverExpressedInteractions(cellchat)

## The number of highly variable ligand-receptor pairs used for signaling inference is 697

# project gene expression data onto PPI (Optional: when running it, USER should set `raw.use = FALSE` in the function `computeCommunProb()` in order to use the projected data)

CellChat infers biologically significant cell-cell communication by assigning each interaction with a probability value and peforming a permutation test. This assigns p-values to the interactions.

By default, CellChat uses type = "triMean" for the computeCommunProbPathway function. This is essentially truncated mean at 25% of the data, which means that the average gene expression is zero if the percent of expressed cells in one group is less than 25%.

However, we can also set the truncated mean value to be something lower. For example, to use 10% truncated mean, we can set type = "truncatedMean" and trim = 0.1.

We can also find out the average expression of different features across cell types using computeAveExpr at different values of truncated mean.

computeAveExpr(cellchat, features = c("CXCL12","CXCR4"), type =  "truncatedMean", trim = 0.25)

##               B       cDC       CTL earlyEry        GMP HSC lateEry Malignant
## CXCR4 0.1661262 0.2331495 0.2612765        0 0.06729584   0       0         0
##             Mono NK       pDC     Plasma      ProB Prog   ProMono          T
## CXCR4 0.07051209  0 0.2132761 0.01283606 0.6419723    0 0.5039044 0.02755188

print("-------------------------------------------------")

## [1] "-------------------------------------------------"

computeAveExpr(cellchat, features = c("CXCL12","CXCR4"), type =  "truncatedMean", trim = 0.1)

##               B      cDC       CTL    earlyEry       GMP        HSC lateEry
## CXCR4 0.3132056 0.333278 0.3892656 0.002062938 0.1772864 0.01140462       0
##       Malignant      Mono        NK       pDC    Plasma      ProB      Prog
## CXCR4         0 0.2006823 0.0941311 0.3319048 0.1417801 0.6653252 0.1285125
##         ProMono         T
## CXCR4 0.5239052 0.1471068

Now, as we can see that most of our cells get zeroed out at a truncated mean of 25%, we will use a truncated mean of 10% for the probability calculation. However, a truncated mean of 25% is more stringent so we get fewer but stronger interactions.

We can filter out the cell-cell communication using filterCommunication if there are only few cells in certain cell groups. By default, the minimum number of cells required in each cell group for cell-cell communication is 10. However, we are going to use a more stringent filter of min.cells = 100.

cellchat = computeCommunProb(cellchat, 
                             type = "truncatedMean",
                             trim = 0.25)

## truncatedMean is used for calculating the average gene expression per cell group. 
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2025-03-18 18:53:37.868048]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2025-03-18 18:55:31.855636]"

#> triMean is used for calculating the average gene expression per cell group. 
#> [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-02-14 00:32:35.767285]"
#> [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2024-02-14 00:33:13.121225]"

cellchat = filterCommunication(cellchat, min.cells = 30)

3.4. Extracting Interaction Information

The function subsetCommunication can be used to access the inferred cell-cell communications of interest. It returns all the inferred cell-cell communications at the level of ligands/receptors. Set slot.name = "netP" to access the the inferred communications at the level of signaling pathways.

interaction_information = subsetCommunication(cellchat)

pathway_information = subsetCommunication(cellchat, 
                                          slot.name = "netP") 

#df.net <- subsetCommunication(cellchat, sources.use = c(1,2), targets.use = c(4,5)) #gives the inferred cell-cell communications sending from cell groups 1 and 2 to cell groups 4 and 5.

#df.net <- subsetCommunication(cellchat, signaling = c("WNT", "TGFb")) #gives the inferred cell-cell communications mediated by signaling WNT and TGFb.

Here we are calculating the communication probability on signaling pathway level by summarizing the communication probabilities of all ligand-receptors interactions associated with each signaling pathway.

cellchat = computeCommunProbPathway(cellchat)

3.5. Aggregating the Network

CellChat calculates the aggregated cell-cell communication network by counting the number of links or summarizing the communication probability. Users can also calculate the aggregated network among a subset of cell groups by setting sources.use and targets.use.

cellchat = aggregateNet(cellchat)

groupSize = as.numeric(table(cellchat@idents)) #This is the number of cells of each type. 

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/interaction_counts.png",
 #   width = 12, height = 12, units = "in", res = 300, bg = "white")

netVisual_circle(cellchat@net$count, 
                 vertex.weight = groupSize, 
                 weight.scale = T, 
                 label.edge= F, 
                 title.name = "Number of interactions")

#dev.off()

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/interaction_weight.png",
 #   width = 12, height = 12, units = "in", res = 300, bg = "white")

netVisual_circle(cellchat@net$weight, 
                 vertex.weight = groupSize, 
                 weight.scale = T, 
                 label.edge= F, 
                 title.name = "Interaction weights/strength")

#dev.off()

Due to the complicated cell-cell communication network, we can examine the signaling sent from each cell group. Here we also control the parameter edge.weight.max so that we can compare edge weights between differet networks.

So, in this step, we are simplifying the network.

Here’s an explanation for what’s going on in the plots:

• The size of the dot on the names of each cell represents the amount of cells of that type.

mat = cellchat@net$weight
#par(mfrow = c(3,4), xpd=TRUE)

for (i in 1:3) 
{
  mat2 = matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] = mat[i, ]
  
  #png(paste0("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/Individual/CellChat_Network_", 
   #          rownames(mat)[i], ".png"), width = 12, height = 12, units = "in", res = 300, bg = "white")

  netVisual_circle(mat2, 
                   vertex.weight = groupSize, 
                   weight.scale = T, 
                   edge.weight.max = max(mat), #The edge size is scaled by the maximum weight.
                   title.name = rownames(mat)[i])
  
  #dev.off()
}

#interaction_information

We can also visualize specific pathways using the signaling parameter of the netVisual_aggregate function.

The graph below shows in the BLTA pathway (from interaction_information$pathway_name). We get the cells that are sources and targets. Arrows go from sources to target (e.g. the ligand BLTA goes from B to T cells and B to NK cells).

We can also access the pathway infromation from cellchat@netP$pathways.

print("Pathways where Malignant cells are the Target: ")

## [1] "Pathways where Malignant cells are the Target: "

unique(interaction_information$pathway_name[interaction_information$target == "Malignant"])

## [1] "MIF"      "RESISTIN" "THBS"     "APP"      "CD99"     "PECAM1"   "SELL"

print("Pathways where Malignant cells are the Source: ")

## [1] "Pathways where Malignant cells are the Source: "

unique(interaction_information$pathway_name[interaction_information$source == "Malignant"])

## [1] "MIF"     "BAFF"    "ANNEXIN" "CypA"    "CD99"    "CLEC"    "MHC-I"  
## [8] "PECAM1"  "PECAM2"

pathways.show = c("APP")

# Hierarchy plot
# Here we define `vertex.receive` so that the left portion of the hierarchy plot shows signaling to fibroblast and the right portion shows signaling to immune cells 
vertex.receiver = seq(1,4) # a numeric vector. 

#png(paste0("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/", pathways.show, ".png"),
 #   width = 12, height = 12, units = "in", res = 300, bg = "white")

netVisual_aggregate(cellchat, 
                    signaling = pathways.show,  
                    vertex.receiver = vertex.receiver)

#dev.off()
# Circle plot
#par(mfrow=c(1,1))
netVisual_aggregate(cellchat,
                    signaling = pathways.show, 
                    layout = "circle")

Here we are visualizing this with a chord figure. The width of the edge represents the strength of the interaction (interaction_information$prob).

#png(paste0("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/", pathways.show, "_chord.png"),
    #width = 12, height = 12, units = "in", res = 300, bg = "white")

netVisual_individual(cellchat, signaling = pathways.show, layout = "chord")

## [[1]]
## 
## [[2]]

#dev.off()

#png(paste0("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/", pathways.show, "_chord_gene.png"),
 #   width = 12, height = 12, units = "in", res = 300, bg = "white")

netVisual_chord_gene(cellchat, 
                     signaling = pathways.show)

## You may try the function `netVisual_chord_cell` for visualizing individual signaling pathway

#dev.off()

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/heatmap.png",
 #   width = 12, height = 8, units = "in", res = 300, bg = "white")

netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")

## Do heatmap based on a single object

#> Do heatmap based on a single object

#dev.off()

3.6. Subsetting the CellChat Object

We can also subset the cellchat to specific cell types using the subsetCellChat function. Here we can pass specific idents using the idents.use parameter.

And then we can plot the chord chart for only these cells.

netVisual_individual(subsetCellChat(cellchat, 
               idents.use = c("B", "T", "Plasma")), 
               signaling = "MIF", 
               layout = "chord")

## The subset of cell groups used for CellChat analysis are  B Plasma T 
## Update slots object@images, object@net, object@netP in a single dataset...

## [[1]]

3.7. Other Analsyes

3.7.1. Contribution Analysis

netAnalysis_contribution(cellchat, signaling = cellchat@netP$pathways[2])

#interaction_information

3.7.2. Specific Bubble Plots

Here we can define bubble plots for all ligand-receptor interactions for cell pairs. The sources and targets are defined according to the levels of idents of the cellchat object (levels(cellchat@idents)). For instance, in our case, we want to understand how GMP, ProMono, and Monocytes are interacting with each other. So, we will classify the sources.use = 5 (factor level of GMP), and the targets.use = c(8, 14).

This gives us all the ligand-receptor pairs, and their interaction strengths (communication probabilities) between GMP and ProMono, and GMP and Monocytes.

levels(cellchat@idents)

##  [1] "B"         "cDC"       "CTL"       "earlyEry"  "GMP"       "HSC"      
##  [7] "lateEry"   "Malignant" "Mono"      "NK"        "pDC"       "Plasma"   
## [13] "ProB"      "Prog"      "ProMono"   "T"

sources = c(1, 2, 3, 10, 16)
targets = c(8)

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/bubble_plot.png",
 #   width = 10, height = 6, units = "in", res = 300, bg = "white")

netVisual_bubble(cellchat, 
                 sources.use = sources, 
                 targets.use = targets,
                 remove.isolate = FALSE)

## Comparing communications on a single object

#dev.off()

Instead of looking at all the pathways, we can specify specific signalling pathways using signaling = c("MIF").

netVisual_bubble(cellchat, 
                 sources.use = sources, 
                 targets.use = targets,
                 signaling = c("PECAM1", "MIF", "APP"),
                 remove.isolate = FALSE)

## Comparing communications on a single object

We can also do the same thing in reverse where we first extract all the enriched receptor-ligand pairs for a specific pathway using extractEnrichedLR and then visualize the enriched pairs using the pairLR.use parameter. The input to this parameter has to be a dataframe with a single column named “interaction_name”.

pairLR.use = extractEnrichedLR(cellchat, signaling = c("MIF"))

netVisual_bubble(cellchat, 
                 sources.use = sources, 
                 targets.use = targets, 
                 pairLR.use = pairLR.use, 
                 remove.isolate = TRUE)

## Comparing communications on a single object

3.7.3. Specific Chord Analysis

We can do the same thing with chord diagrams where we can specify sources using sources.use and targets using targets.use (based on the levels of idents of the cellchat object - levels(cellchat@idents)).

And here we can also specify specific signaling pathways (cellchat@netP$pathways)to look at like MIF for instance.

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/CD99_Chord.png",
 #   width = 7, height = 7, units = "in", res = 300, bg = "white")

netVisual_chord_cell(cellchat, 
                     sources.use = sources, 
                     targets.use = targets,
                     signaling = c("CD99"),
                     lab.cex = 0.5,
                     legend.pos.y = 30)

## Plot the aggregated cell-cell communication network at the signaling pathway level

#dev.off()

Or we can specify multiple sources and a single target.

netVisual_chord_cell(cellchat, 
                     sources.use = sources,
                     targets.use = targets,
                     signaling = "CD99",
                     lab.cex = 0.5,
                     legend.pos.y = 30)

## Plot the aggregated cell-cell communication network at the signaling pathway level

3.7.4. Violin Plots

We can also plot violin plots of the expression levels of all ligands (interaction_information$ligand) and receptors (interaction_information$receptor) that are enriched in a specific pathway (interaction_information$pathway_name) across all cell types.

unique(interaction_information$pathway_name)

##  [1] "MIF"      "BAFF"     "RESISTIN" "ANNEXIN"  "CTSG"     "CypA"    
##  [7] "THBS"     "APP"      "CD99"     "CLEC"     "MHC-I"    "PECAM1"  
## [13] "SELL"     "PECAM2"

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/violin_plot.png",
 #   width = 18, height = 12, units = "in", res = 300, bg = "white")

plotGeneExpression(cellchat, 
                   signaling = c("MIF", "PECAM1"), 
                   enriched.only = TRUE, 
                   type = "violin")

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

#dev.off()

3.7.5. Dominant Senders and Receivers

cellchat = netAnalysis_computeCentrality(cellchat, slot.name = "netP") # the slot 'netP' means the inferred intercellular communication network of signaling pathways

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/sender_strength_plot.png",
 #   width = 9, height = 4, units = "in", res = 300, bg = "white")

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(cellchat, 
                                  signaling = c("MIF"),
                                  width = 16, 
                                  height = 4, 
                                  font.size = 10)

#dev.off()

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/incoming_outgoing_strength_plot.png",
 #   width = 15, height = 7, units = "in", res = 300, bg = "white")

# Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
gg1 = netAnalysis_signalingRole_scatter(cellchat)

## Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

# Signaling role analysis on the cell-cell communication networks of interest
gg2 = netAnalysis_signalingRole_scatter(cellchat, signaling = c("MIF"))

## Signaling role analysis on the cell-cell communication network from user's input

gg1 + gg2

#dev.off()

Here we can see that for MIF pathway, the GMP cells are the dominant senders compared to ProMono, that are the dominant receivers.

gg2

We can also use heatmaps to plot the interaction strenth of individual ligands and cell types for different signaling pathways.

The top colored bar plot shows the total signaling strength of a cell group by summarizing all signaling pathways displayed in the heatmap. The right grey bar plot shows the total signaling strength of a signaling pathway by summarizing all cell groups displayed in the heatmap.

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/incoming_outgoing_heatmap.png",
 #   width = 12, height = 8, units = "in", res = 300, bg = "white")

ht1 = netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", height = 15, width = 10)
ht2 = netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", height = 15, width = 10)
ht1 + ht2

#dev.off()

And of course, like everything else, we can do it for a specific pathway.

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/incoming_outgoing_pathway_heatmap.png",
 #   width = 7, height = 6, units = "in", res = 300, bg = "white")

ht = netAnalysis_signalingRole_heatmap(cellchat, 
                                       signaling = c("CD99", "MIF"),
                                       pattern = "outgoing",
                                       width = 14,
                                       height = 8)
ht

#dev.off()

3.7.6. Signaling Coordination Analysis

We can also see how the different incoming and outgoing singals are coordinated by the cells using the identifyCommunicationPatterns function. The correct paramter k can be selected using the selectK function, however, it takes a very long time to run.

How to choose the best k: selectK(cellchat, pattern = "outgoing")

#selectK(cellchat, pattern = "outgoing")

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/pattern_heatmap.png",
 #   width = 14, height = 12, units = "in", res = 300, bg = "white")

cellchat = identifyCommunicationPatterns(cellchat, 
                                         pattern = "outgoing", 
                                         k = 2,
                                         height = 17,
                                         width = 9)

#dev.off()

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/pattern_flow.png",
 #   width = 14, height = 12, units = "in", res = 300, bg = "white")

netAnalysis_river(cellchat, pattern = "outgoing")

## Please make sure you have load `library(ggalluvial)` when running this function

#dev.off()

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/pattern_bubble.png",
 #   width = 14, height = 8, units = "in", res = 300, bg = "white")

netAnalysis_dot(cellchat, pattern = "outgoing")

#dev.off()

3.7.7. Signaling Similarity Analysis

Further, CellChat is able to quantify the similarity between all significant signaling pathways and then group them based on their cellular communication network similarity. Grouping can be done either based on the functional or structural similarity.

Functional similarity: High degree of functional similarity indicates major senders and receivers are similar, and it can be interpreted as the two signaling pathways or two ligand-receptor pairs exhibit similar and/or redundant roles. The functional similarity analysis requires the same cell population composition between two datasets.

Structural similarity: A structural similarity was used to compare their signaling network structure, without considering the similarity of senders and receivers.

#cellchat = readRDS("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/cellchat.rds")

set.seed(1010)

cellchat = computeNetSimilarity(cellchat, 
                                type = "functional")

cellchat = netEmbedding(cellchat, 
                        type = "functional")

## Manifold learning of the signaling networks for a single dataset

cellchat = netClustering(cellchat, 
                         type = "functional",
                         methods = "kmeans",
                         k = 2)

## Classification learning of the signaling networks for a single dataset

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/umap_functional.png",
 #   width = 10, height = 7, units = "in", res = 300, bg = "white")

# Visualization in 2D-space
netVisual_embedding(cellchat, 
                    type = "functional", 
                    label.size = 3.5)

#dev.off()

We can also plot the groups individually

netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)

set.seed(1010)

cellchat = computeNetSimilarity(cellchat, 
                                type = "structural")

cellchat = netEmbedding(cellchat, 
                        type = "structural")

## Manifold learning of the signaling networks for a single dataset

cellchat = netClustering(cellchat, 
                         type = "structural",
                         methods = "kmeans",
                         k = 2)

## Classification learning of the signaling networks for a single dataset

#png("/Users/ali/Desktop/NYU/4th Semester/Single-cell Omics/Second Homework/umap_structural.png",
 #   width = 10, height = 7, units = "in", res = 300, bg = "white")

# Visualization in 2D-space
netVisual_embedding(cellchat, 
                    type = "structural", 
                    label.size = 3.5)

#dev.off()